get---
title: "R_SM_A5"
author: "Loreto Cox & Carmen Le Foulon"
date: "2024-01-26"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(brms)
library(tidyverse)
library(tidybayes)

```

# Info

Replication: Supplementary Appendix 5

We follow Tim Mastny’s code available at https://timmastny.com/blog/multilevel-mrp-tidybayes-brms-stan/.

# Data


## CENSUS

Created datasets from 2017 Census Data from Chilean Instituto Nacional de Estadísticas (INE)

In 2018, some counties changed their district, we created two different population data bases, one for  2017 and before and the other, for the 2021 election. 

See "Codebook.docx" for more details on the two datasets used: "pop_elec2017.rds" and "pop_elec2021.rds"

## District magnitude

Dataset with information on post-reform district id and district magnitude: "m_dip.rds", see "Codebook.docx" for more details.


## Centro de Estudios Público (CEP) Public Opinion Surveys

In order to assign the district to each respondent in each survey, we used municipality data which is not publicly available. Hence, we have made available the reduced dataset with information on the recoded key variables used in the analysis and the district variable.

The data (except municipality variable) is the publicly available in the consolidated dataset avaiable at https://www.cepchile.cl/opinion-publica/encuesta-cep/ downloaded on June 15, 2023.


# Estimation
## CEP 69

```{r}
ncep<- 69
nyear <- 2013
nelec <- 2013

pop <- readRDS("data/SM_5/cep_data/pop_elec2017.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep69mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_69",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_69.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_69",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_69",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_69.rds")

```


### Democracy is working bad o very bad


```{r}
m_evaldem_na <- brm( evaldem_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_evaldem_na_69",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_evaldem_na)


estim_evaldem_na <- 
  m_evaldem_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_evaldem_na <- 
  estim_evaldem_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "evaldem_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_evaldem_na <- estim_dist_evaldem_na %>% select(-m_dip, -tipo)
colnames(estim_dist_evaldem_na) <- nomcols

 tot_evaldem <-  estim_dist_evaldem_na %>% 
        mutate(cep=ncep, elec=nelec, year=nyear) 

saveRDS(tot_evaldem , "data/SM_5/output_results/estim_evaldem_69.rds")
```



## CEP 71

```{r}
ncep<- 71
nyear <- 2014
nelec <- 2013


pop <- readRDS("data/SM_5/cep_data/pop_elec2017.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep71mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))
```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_71",
                    control = list(adapt_delta = 0.98),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_71.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_71",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_71",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_71.rds")

```



## CEP 79

```{r}
ncep<- 79
nyear <- 2017
nelec <- 2017

pop <- readRDS("data/SM_5/cep_data/pop_elec2017.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep79mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))


```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_79",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_79.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_79",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_79",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_79.rds")

```


### Democracy is working bad o very bad


```{r}
m_evaldem_na <- brm( evaldem_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_evaldem_na_79",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_evaldem_na)


estim_evaldem_na <- 
  m_evaldem_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_evaldem_na <- 
  estim_evaldem_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "evaldem_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_evaldem_na <- estim_dist_evaldem_na %>% select(-m_dip, -tipo)
colnames(estim_dist_evaldem_na) <- nomcols

 tot_evaldem <-  estim_dist_evaldem_na %>% 
        mutate(cep=ncep, elec=nelec, year=nyear) 

saveRDS(tot_evaldem , "data/SM_5/output_results/estim_evaldem_79.rds")
```



## CEP 82

```{r}
ncep<- 82
nyear <- 2018
nelec <- 2017


pop <- readRDS("data/SM_5/cep_data/pop_elec2017.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep82mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))
```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_82",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_82.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_82",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_82",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_82.rds")

```




## CEP 85

```{r}
ncep<- 85
nyear <- 2021
nelec <- 2021

pop <- readRDS("data/SM_5/cep_data/pop_elec2021.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep85mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))


```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_85",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_85.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_85",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_85",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_85.rds")

```


### Democracy is working bad o very bad


```{r}
m_evaldem_na <- brm( evaldem_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_evaldem_na_85",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_evaldem_na)


estim_evaldem_na <- 
  m_evaldem_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_evaldem_na <- 
  estim_evaldem_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "evaldem_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_evaldem_na <- estim_dist_evaldem_na %>% select(-m_dip, -tipo)
colnames(estim_dist_evaldem_na) <- nomcols

 tot_evaldem <-  estim_dist_evaldem_na %>% 
        mutate(cep=ncep, elec=nelec, year=nyear) 

saveRDS(tot_evaldem , "data/SM_5/output_results/estim_evaldem_85.rds")
```



## CEP 86

```{r}
ncep<- 86
nyear <- 2022
nelec <- 2021

pop <- readRDS("data/SM_5/cep_data/pop_elec2021.rds")
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds") %>% group_by(m_dip) %>% mutate(tipo=row_number()) %>% 
            ungroup() %>% mutate(distrito=factor(distrito))
cep<- readRDS("data/SM_5/cep_data/cep86mrp_dist.rds")
```


- Create factor variables
```{r}
cep <- cep %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))

pop <- pop %>%  mutate(educ=factor(educ), edad=factor(edad), urb=factor(urb),
                      distrito=factor(distrito))
```

We follow recommendations from Stan Development Team(Last updated: 2022-03-10) available at:
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

### Distrust in Congress

```{r}


m_noconf_cong_na <- brm( noconf_cong_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                    file = "data/SM_5/output_mrp/m_noconf_cong_na_86",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_noconf_cong_na)

estim_noconf_cong_na <- 
  m_noconf_cong_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_noconf_cong_na <- 
  estim_noconf_cong_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()

nombre <- "noconf_cong_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_noconf_cong_na <- estim_dist_noconf_cong_na %>% select(-m_dip, -tipo)
colnames(estim_dist_noconf_cong_na) <- nomcols

 tot_conf_cong <- estim_dist_noconf_cong_na %>% 
                  mutate(cep=ncep, elec=nelec, year=nyear) 

 saveRDS(tot_conf_cong, "data/SM_5/output_results/estim_conf_cong_86.rds")

```

### Country is in decadence & bad or very bad conomic conditions 

```{r}

## Country is in decadence
m_endecad_na <- brm( endecad_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_endecad_na_86",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_endecad_na)

estim_endecad_na <- 
  m_endecad_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_endecad_na <- 
  estim_endecad_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


## Economic conditions of the country are bad o very bad


m_econmal_na <- brm( econmal_na ~ muj + (1|urb) + (1|edad) + (1|educ) + (1|distrito), 
                    data = cep, 
                    family = bernoulli(),
                   file = "data/SM_5/output_mrp/m_econmal_na_86",
                    control = list(adapt_delta = 0.95),
                    prior=c(set_prior("normal(0,0.2)", class='b'),
                          set_prior("normal(0,0.2)", class='sd', group="urb"),
                            set_prior("normal(0,0.2)", class='sd', group="educ"),
                            set_prior("normal(0,0.2)", class='sd', group="edad"),
                            set_prior("normal(0,0.2)", class='sd', group="distrito")))

summary(m_econmal_na)


estim_econmal_na <- 
  m_econmal_na %>% 
  tidybayes::add_predicted_draws(newdata = pop) %>% 
  rename(alp_predict = .prediction) %>% 
  mutate(alp_predict_prop = alp_predict*per_dist) %>% 
  group_by(distrito, .draw) %>% 
  summarise(alp_predict = sum(alp_predict_prop)) %>%
  ungroup()

estim_dist_econmal_na <- 
  estim_econmal_na %>% 
  group_by(distrito) %>% 
  summarise(mean = mean(alp_predict), 
            lower = quantile(alp_predict, 0.025), 
            upper = quantile(alp_predict, 0.975)) %>% 
  ungroup() %>% full_join(mdip) %>% 
  group_by(m_dip) %>% mutate(tipo=row_number()) %>% ungroup()


nombre <- "endecad_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_endecad_na <- estim_dist_endecad_na %>% select(-m_dip, -tipo)
colnames(estim_dist_endecad_na) <- nomcols

nombre <- "econmal_na"
nomcols <-c("distrito", paste(c("mean", "lower", "upper"), nombre, sep="_"))
estim_dist_econmal_na <- estim_dist_econmal_na %>% select(-m_dip, -tipo)
colnames(estim_dist_econmal_na) <- nomcols


 tot_mb <- full_join(estim_dist_endecad_na, estim_dist_econmal_na) %>%
            mutate(cep=ncep, elec=nelec, year=nyear) 
 
 
saveRDS(tot_mb, "data/SM_5/output_results/estim_mb_86.rds")

```

# Figures

## Figure SM4: Estimates for distrust in Congress

```{r}

mdip <- readRDS("data/SM_5/cep_data/m_dip.rds")

c69<-  readRDS("data/SM_5/output_results/estim_conf_cong_69.rds")
c71<-  readRDS("data/SM_5/output_results/estim_conf_cong_71.rds")
c79<-  readRDS("data/SM_5/output_results/estim_conf_cong_79.rds")
c82<-  readRDS("data/SM_5/output_results/estim_conf_cong_82.rds")
c85<-  readRDS("data/SM_5/output_results/estim_conf_cong_85.rds")
c86<-  readRDS("data/SM_5/output_results/estim_conf_cong_86.rds")

conf <- bind_rows(c69, c71, c79, c82, c85, c86) %>% full_join(mdip)


conf <- conf %>% mutate(pre=ifelse(cep==69 | cep==79 | cep==85, "Pre-election",
                          "Post-election"),
                        pre=factor(pre, ordered=T, levels=c("Pre-election", "Post-election")))

ggplot(conf) + geom_smooth(aes(y=mean_noconf_cong_na, 
                             x=m_dip, 
                             colour=pre), se=FALSE, linewidth=1.0)   +
  geom_jitter(aes(y=mean_noconf_cong_na, x=m_dip, colour=pre,
                  shape=pre), 
              width = 0.08, size=2) +
  facet_grid(cols =vars(elec)) +
  theme_bw() + scale_y_continuous(limits=c(.5,1), expand = c(0,0))  +
  
  scale_colour_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( "darkgoldenrod1", "darkorchid4")) +
  scale_shape_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( 19, 17)) +
    labs(caption = "Data source: CEP Surveys 69, 71, 79, 82, 85, 86 ") +

  xlab("District Magnitude") + ylab("Percent") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size = 10),
        title =element_text(size=10),
        legend.text=element_text(size=10),
        legend.position = "bottom")
```

## Figure SM5: Estimates for perception that the country is in decadence

```{r}
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds")
c69<-  readRDS("data/SM_5/output_results/estim_mb_69.rds")
c71<-  readRDS("data/SM_5/output_results/estim_mb_71.rds")
c79<-  readRDS("data/SM_5/output_results/estim_mb_79.rds")
c82<-  readRDS("data/SM_5/output_results/estim_mb_82.rds")
c85<-  readRDS("data/SM_5/output_results/estim_mb_85.rds")
c86<-  readRDS("data/SM_5/output_results/estim_mb_86.rds")

mb <- bind_rows(c69, c71, c79, c82, c85, c86) %>% full_join(mdip)

mb <- mb %>% mutate(pre=ifelse(cep==69 | cep==79 | cep==85,
                               "Pre-election","Post-election"),
                        pre=factor(pre, ordered=T,
                        levels=c("Pre-election", "Post-election")))

ggplot(mb) + geom_smooth(aes(y=mean_endecad_na, 
                             x=m_dip, 
                             colour=pre),se=FALSE, linewidth=1.0)   +
  geom_jitter(aes(y=mean_endecad_na, x=m_dip, colour=pre,
                  shape=pre), 
              width = 0.08, size=2) +
  facet_grid(cols =vars(elec)) +
  theme_bw() + scale_y_continuous(limits=c(0,.5), expand = c(0,0))  +
  
  scale_colour_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( "darkgoldenrod1", "darkorchid4")) +
  scale_shape_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( 19, 17)) +
    labs(caption = "Data source: CEP Surveys 69, 71, 79, 82, 85, 86 ") +

  xlab("District Magnitude") + ylab("Percent") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size = 10),
        title =element_text(size=10),
        legend.text=element_text(size=10),
        legend.position = "bottom")
```

## Figure SM6: Estimates for perception that economic conditions of the country are bad o very bad

```{r}
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds")
c69<-  readRDS("data/SM_5/output_results/estim_mb_69.rds")
c71<-  readRDS("data/SM_5/output_results/estim_mb_71.rds")
c79<-  readRDS("data/SM_5/output_results/estim_mb_79.rds")
c82<-  readRDS("data/SM_5/output_results/estim_mb_82.rds")
c85<-  readRDS("data/SM_5/output_results/estim_mb_85.rds")
c86<-  readRDS("data/SM_5/output_results/estim_mb_86.rds")

mb <- bind_rows(c69, c71, c79, c82, c85, c86) %>% full_join(mdip)

mb <- mb %>% mutate(pre=ifelse(cep==69 | cep==79 | cep==85,
                               "Pre-election","Post-election"),
                        pre=factor(pre, ordered=T,
                        levels=c("Pre-election", "Post-election")))

ggplot(mb) + geom_smooth(aes(y=mean_econmal_na, 
                             x=m_dip, 
                             colour=pre), se=FALSE, linewidth=1.0)   +
  geom_jitter(aes(y=mean_econmal_na, x=m_dip, colour=pre,
                  shape=pre), 
              width = 0.08, size=2) +
  facet_grid(cols =vars(elec)) +
  theme_bw() + scale_y_continuous(limits=c(0,.8), expand = c(0,0))  +
  
  scale_colour_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( "darkgoldenrod1", "darkorchid4")) +
  scale_shape_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( 19, 17)) +
    labs(caption = "Data source: CEP Surveys 69, 71, 79, 82, 85, 86 ") +

  xlab("District Magnitude") + ylab("Percent") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size = 10),
        title =element_text(size=10),
        legend.text=element_text(size=10),
        legend.position = "bottom")
```

## Figure SM7: Estimates for perception that democracy is working bad o very bad

```{r}
mdip <- readRDS("data/SM_5/cep_data/m_dip.rds")
c69<-  readRDS("data/SM_5/output_results/estim_evaldem_69.rds")
c79<-  readRDS("data/SM_5/output_results/estim_evaldem_79.rds")
c85<-  readRDS("data/SM_5/output_results/estim_evaldem_85.rds")

eval <- bind_rows(c69,  c79,  c85) %>% full_join(mdip)

eval <- eval %>% mutate(pre=ifelse(cep==69 | cep==79 | cep==85,
                               "Pre-election","Post-election"),
                        pre=factor(pre, ordered=T,
                        levels=c("Pre-election", "Post-election")))


ggplot(eval) + geom_smooth(aes(y=mean_evaldem_na, 
                             x=m_dip, 
                             colour=pre), se=FALSE, linewidth=1.0)   +
  geom_jitter(aes(y=mean_evaldem_na, x=m_dip, colour=pre,
                  shape=pre), 
              width = 0.08, size=2) +
  facet_grid(cols =vars(elec)) +
  theme_bw() + scale_y_continuous(limits=c(0,.5), expand = c(0,0))  +
  
  scale_colour_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( "darkgoldenrod1", "darkorchid4")) +
  scale_shape_manual(name = "Survey", labels = c("Pre-election", "Post-election"), values = c( 19, 17)) +
  
  xlab("District Magnitude") + ylab("Percent") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size = 10),
        title =element_text(size=10),
        legend.text=element_text(size=10),
        legend.position = "bottom")
```

